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ABSTRACT 


A  modal  approximation  for  the  self  and  mutual  radiation 
impedances  has  been  derived  for  arrays  of  spherical 
transducers  with  small  ka  values,  where  ka  is  the  acoustic 
wave . number  multiplied  by  the  radius  of  the  sphere .  I  term 
this  the  "Modal  Pritchard"  approximation,  as  it  is  related 
to  the  so-called  Pritchard  approximation,  often  employed  to 
calculate  mutual  radiation  impedances .  I  investigated  the 
utility  of  the  approximate  mutual  radiation  impedance 
expression  for  three  two-body  problems  (monopoles,  aligned 
dipoles,  and  aligned-linear  quadrupoles) .  For  these  cases, 
approximate  values  were  found  to  be  in  good  agreement  with 
those  obtained  using  a  full  Spherical  Addition  Theorem 
calculation,  and  are  an  improvement  over  the  simple 
Pritchard  approximation.  Additionally,  I  investigated  the 
mutual  radiation  impedance  expression  in  one  particular 
three-body  problem.  Because  of  the  Modal  Pritchard 
approximation's  inability  to  correctly  handle  scattering,  we 
recommend  using  the  full  Spherical  Addition  Theorem 
calculation  when  scattering  is  important. 

Finally,  I  investigated  the  use  of  a  new  finite  element 
mesh  to  calculate  the  T-matrix  for  a  given  transducer.  The 
T-matrix  relates  the  incident  and  scattered  waves  for  a 
single  transducer,  in  an  orthogonal  (spherical  harmonic) 
basis  set.  The  monopole  element  showed  an  increase  in 
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error,  while  we  saw  some  improvements  in  the  higher-order 
diagonal  elements.  Off-diagonal  elements,  which  should  be 
zero  for  a  spherical  scatter,  were  satisfactorily  small  in 
most  cases.  Although  the  results  were  less  than  favorable, 
I  was  able  to  streamline  the  T-matrix  calculation  while 
providing  a  new  method  of  examining  the  off-diagonal 
elements . 
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I.  INTRODUCTION 

The  work  described  in  this  thesis  is  part  of  an  ongoing 
effort  to  improve  our  ability  to  accurately  model  densely 
packed  underwater  acoustic  arrays.  When  an  acoustic  array's 
operating  characteristics  are  closely  predicted  we  say  the 
array  behaves  well.  This  behavior  is  influenced  by  element 
interaction,  characterized  by  a  parameter  called  mutual 
radiation  impedance  (Ref.  1).  A  worst  case  example  of  "bad" 
behavior  is  that  a  transducer  element  may  actually  absorb 
acoustic  power,  which  is  then  absorbed  back  into  the  driving 
amplifier,  thus  possibly  causing  an  overload  and  electrical 
failure  (Ref.  2)  .  Other  examples  of  bad  behavior  include 
changes  in  the  radiated  source  level  and  deviations  of  the 
beampattern,  steering  direction,  and  side  lobe  amplitudes. 

Element- to-element  interactions  are  more  significant  in 
arrays  that  are  packed  into  a  volume  of  small  size,  which  is 
the  trend  in  low  frequency  active  arrays.  We  will  analyze 
the  mutual  radiation  impedance  for  spherical  sources  of 
small  ka  values,  where  k  is  the  wave  number  and  a  is  the 
radius  of  the  sphere.  The  mutual  radiation  impedance  is  an 
important  parameter  used  to  calculate  the  pressure  on  one 
array  transducer  due  to  the  radiation  and  scattering  from 
other  transducers  in  the  array. 

Additionally,  I  will  use  a  finite  element  code  to  solve 
for  the  radiating  and  scattering  properties  of  a  single 
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transducer.  Therefore,  this  thesis  consists  of  two 
portions . 

The  first  part  is  devoted  to  deriving  a  general 
expression  for  the  mutual  radiation  impedance  between  two 
identical  radiating  hard  spheres  in  an  otherwise  free  field 
environment  using  a  normal  mode  approach.  Additionally,  an 
expression  for  the  self-radiation  impedance,  which  is  the 
radiation  loading  on  a  single  element  in  a  free  field,  is 
found. 

I  have  compared  our  modal  approximation  for  the  mutual 
radiation  impedance  with  the  commonly-used  so-called 
"Pritchard"  approximation  for  the  particular  case  ka  =  1, 
and  will  show  the  improvement  possible  using  our  modal 
approximation.  This  improvement  is  particularly  noticeable 
as  the  center-to-center  distance  between  sources  is  reduced. 
In  this  "Simple  Pritchard"  approximation,  as  I  refer  to  it, 
the  real  part  of  the  self-radiation  impedance  is  multiplied 
by  a  spherical  Hankel  function  of  zeroth  order  to  find  the 
mutual  radiation  impedance .  I  will  show  how  this  "simple" 
approximation  starts  to  break  down  as  the  transducers  are 
spatially  brought  closer  together. 

The  Addition  Theorem  (Ref.  3)  provides  one  the  ability 
to  mathematically  express  the  pressure  radiated  from  one 
body  referenced  to  the  coordinate  system  of  another  body. 
Using  the  Addition  Theorem,  it  is  possible  to  very 
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accurately  compute  the  mutual  radiation  impedance  of  a 
sphere  in  an  arbitrarily  densely  packed  array  of  any  number 
of  elements.  I  have  used  the  Addition  Theorem  as  a 
computational  standard. 

I  also  applied  our  modal  approximation  to  a  specified 
three-body  problem  to  further  investigate  our  modal 
approximation,  which  I  will  hereafter  refer  to  as  the  "Modal 
Pritchard"  approximation.  In  this  problem  we  note  the 
inability  of  our  Modal  Pritchard  approximation  to  handle  the 
effects  of  scattering.  Because  of  this  breakdown  we 
recommend  using  the  full  Spherical  Addition  Theorem  for  this 
problem  and  other  cases  when  scattering  becomes  important. 
Methods  and  computer  codes  for  this  portion  of  the  thesis 
are  provided  in  Appendixes  A  and  B. 

The  second  part  of  this  thesis  makes  extensive  use  of  a 
finite  element  code,  named  ATILA,  which  was  written  by 
engineers  at  the  Institut  Superieur  d' Electronique  du  Nord 
(ISEN) .  ATILA  has  been  specifically  developed  to  aid  in  the 
design  of  sonar  transducers,  but  is  also  used  to  model  a 
variety  of  acoustic  problems. 

The  final  portion  of  this  thesis  addresses  the  T-matrix 
method  (Ref.  3) .  In  order  to  model  certain  aspects  of  array 
operation  such  as  near-field  pressure,  it  is  necessary  to 
solve  for  the  radiating  and  scattering  properties  of  a 
single  transducer.  We  accomplish  this  calculation  by  using 
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the  ATILA  finite  element  code.  The  T-Matrix  (or  Transition 
Matrix)  characterizes  the  scattering  properties  of  a 
specific  body  by  translating  incoming  pressure  waves  to 
outgoing  ones. 

I  will  provide  a  brief  description  of  ATILA 's  method  of 
calculation  that  distinguishes  it  from  other  finite  element 
codes.  Readers  less  familiar  with  finite  element  codes  may 
want  to  seek  further  information  from  reference  books  such 
as  Kwon's  "The  Finite  Element  Method  using  MATLAB"  (Ref.  4). 

Previous  work  with  ATILA  to  calculate  the  T-matrix  for 
a  thin- spherical  shell  has  been  shown  to  be  in  error  when 
compared  with  a  verified  analytical  solution  (Ref.  5)  . 
Various  limitations  of  the  ATILA  code  may  explain  the  errors 
present.  I  investigated  whether  improvements  are  obtained 
by  making  some  corrections  to  ATILA' s  code  and  by  the 
implementation  of  a  new,  refined,  finite  element  mesh.  This 
new  mesh  increased  the  errors  to  the  monopole  element  but 
did  provide  some  improvement  to  other  T-matrix  elements. 
Further  improvements  may  be  realized  with  the  advent  of  more 
sophisticated  radiation  boundary  operators.  An  early 
delivery  of  an  addition  to  the  ATILA  code,  EQI,  which  uses 
using  a  boundary  element  method  to  handle  the  fluid,  did  not 
occur  in  time  for  use  with  this  thesis.  Investigation  of 
the  improvement  obtained  using  ATILA  coupled  with  EQI  will 
have  to  be  the  subject  of  future  research. 
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In  addition  to  providing  some  improvement  in  the 
calculation  of  the  T-matrix,  I  have  been  able  to  streamline 
the  procedure .  Output  scattered  pressures  from  ATILA  were 
post-processed  using  easy  to  read  MATLAB  program  scripts. 
The  resulting  T-matrix  calculations  were  quickly  compared 
with  the  analytical  results  (also  done  in  MATLAB)  .  Listed 
in  Appendix  C  are  the  MATLAB  codes  used  to  calculate  the 
theoretical  values  of  the  T-matrix  elements  and  the  T-matrix 
elements  calculated  from  ATILA' s  scattered  pressure  results. 

In  summary,  by  deriving  an  improved  analytical 
formulation  for  the  mutual  radiation  impedance,  one  can 
quickly  investigate  the  effects  of  parameter  variations, 
such  as  the  distance  between  elements  in  the  array. 
Additionally  since  this  formulation  involves  angular 
dependence,  I  have  devised  a  means  of  investigating  the 
effects  of  altering  the  orientation  between  elements  in  an 
array.  We  recommend  using  the  Modal  Pritchard  approximation 
for  approximate  solutions  when  scattering  is  unimportant, 
and  the  full  spherical  Addition  Theorem  otherwise. 
Streamlining  the  T-matrix  calculation  for  a  given  transducer 
will  allow  follow-on  work  to  be  accomplished  with  less 
startup  time. 
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II.  THEORETICAL  BACKGROUND  SECTION 

This  chapter  addresses  the  development  of  a  modal 
approximation  for  the  self  and  mutual  radiation  impedances 
for  two  spheres  in  an  otherwise  free  field  environment. 


A.  BOUNDARY  CONDITIONS 

Consider  a  space  and  time  dependent  velocity  potential, 
defined  as  iP"(r,  t)  =  *F(r)  e5<at ,  where  the  vector  r  =  ( r,6 ,(f)) 

represents  the  spatial  dependence. .  A  relationship  between 
the  pressure  (p)  and  the  velocity  potential  ( ¥0  states  p  = 
-po9'P/9t  (Ref.  6),  so  that  p  =  -j(opoW. 

Taking  the  gradient  and  then  the  dot  product  with  the 
radial  unit  vector,  we  obtain  dW/dr  -  -1/ (jcopo)  dp/dr.  By 
noting  that  the  velocity  vector  is  defined  to  equal  the 
gradient  of  the  velocity  potential  (u  =  VW)  ,  we  can  derive 
a  boundary  condition  for  a  surface  normal  velocity.  For  a 
sphere  of  radius  a,  for  example,  evaluation  of  dW/dr  =  - 
1/ ( j(opo)  dp! dr  at  r=a  provides  an  expression  for  the  radial 
velocity  on  the  sphere's  surface.  This  will  be  a  boundary 
condition  for  our  analysis. 


j  fy 
<°Po  ^  r=a 
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B.  PROBLEM  SETUP  AMD  ADDITION  THEOREM 

Consider  two  identical  spheres  (radii  -a)  spaced  with 

a  center-to-center  distance  =  d.  Letting  (rx,  61,  0X)  and 
(r2,02,4>2)  be  the  spherical  coordinates  relative  to  the 
centers  of  spheres  one  and  two,  respectively.  With  the  two 
spheres  radiating  in  distinct  normal  modes,  the  radial 
velocities  are  expressed  below. 

(2) 

Where  £2”  =  p”(cos0)eJm<|,is  a  modal  eigenfunction,  p™  is  the 
Legendre  function,  and  VTiandVn2  are  modal  amplitudes . 

The  pressure  from  each  sphere  can  be  mathematically 
expressed  as  shown  in  Equation  (3)  (Ref.  3)  .  Since  we  are 
assuming  a  eicot  time  dependence,  we  use  a  spherical  Hankel 
function  of  the  second  kind.  If  we  were  to  assume  a  time 
dependence  of  e'iQ* ,  then  we  would  instead  use  a  spherical 
Hankel  function  of  the  first  kind.  After  Equation  (3),  we 
drop  the  Hankel  function's  (2)  superscript. 


Given  the  Addition  Theorem,  we  can  express  the  pressure 
from  sphere  two  relative  to  sphere  one  (see  Equation  (4), 
Ref.  7).  The  appendix  within  Ref.  7  provides  the  definition 
for  the  function  a(v,p,  t,ju,  s)  .  The  angles  012  and  $12  are 
formed  by  describing  the  location  of  sphere  one's  center 
relative  to  sphere  two's  center  (see  Figure  1)  .  The  prime 
in  the  last  summation  means  the  index  is  incremented  by  two. 

/ 

1  i  'f  a(v,.  j v  (kr,)* 

v,=0  ■Mlm-vlpAt1-v ,1 

P?\s2~M  J 

hPl{kd  'P?iel,</>)Q.p-'‘i6n,<l>n) 

1  Vi  *  1 

A  reapplication  of  the  Addition  Theorem  provides  us  a 
mathematical  expression  for  the  pressure  from  sphere  one 

relative  to  sphere  2  coordinates . 

/ 

X  X  X  ^(v2> Pi’ti’Mz’Si)  jv  (kri)* 

p2=l  i  a'  r° 

ri  Lj  Lj  PiVrMzl 

tr°  Sr-U 

hpjjcd  pf'fovfo) 

**  V2  *  2 
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Figure  1.  Orientation  of  sphere's  origins 


C.  APPLICATION  OF  BOUNDARY  CONDITIONS 

Recalling  the  boundary  condition  for  the  sphere's 
surface  velocity  as  stated  in  Equation  (1)  we  obtain 


After  the  application  of  the  boundary  conditions,  we 
use  normal  mode  analysis  (orthogonality)  to  eliminate  some 
terms.  In  Equation  (3)  we  note  that  the  non- zero  terms  are 
those  for  which  t2  =  n2,  s±  =  2T\,  t2  =  n2,  and  s2  =  ir^. 

Similarly,  in  Equations  (4)  and  (5)  the  non- zero  terms  are 
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those  where  v1  -  n2,  ju2  =  it, v2  =  n2,  and  JU2  =  ir^.  Then  by 
using  the  relationship  with  sound  speed,  angular  frequency, 
and  wave  number  (c  =  co/k)  and  dividing  by  a  common  term  we 
obtain 


ip0cVm' 

— - —  =  A1  S  S  + 


h'nm 


~  A  '  f  (ka)  tsnt 
t2~o  s2  —t2  (ka)  p-urn,] 


-ipcVm‘ 

—  *  =  A2  s  5  + 

h'n  ( ka)  nrm  mi  mmi 

00  X?'  7  (^<2)  /  V 

X  X  Ats  ~~f  TT  7  X  (2  \fl2’  P2>  fulfil’ Si/' 

tr 0  *~-f,  llS'h  n2{ka)  p  =  tl-n2\ 

pAsrmh 

h„  (j kd  )or m(x-eaji+0j 


Note  that  the  first  summation  in  each  of  these  last  two 
equations  involves  infinity.  For  computational  purposes  we 
must  truncate  to  some  finite  number  ( K )  .  Note  that  when  K=  0 
the  sphere  is  vibrating  in  the  monopole  mode.  For  K=1 
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dipole  modes  are  kept,  and  when  K=  2,  quadrupole  modes  are 
present,  and  so  on  to  higher  modes. 

D.  SOLVING  FOR  THE  AMPLITUDE  COEFFICIENTS 

We  can  write  these  last  two  equations  in  a  matrix  form 
by  first  constructing  a  large  column  vector  whose  components 
are  the  stacked  amplitude  coefficients  for  spheres  one  and 
two.  The  ordering  of  the  coefficients  for  the  spheres  are 
provided  below,  where  A3=  [A00  AU1  A^  A^  ...  .  ]T  is  the  vector  of 
coefficients  for  sphere  one. 


The  C  vectors  above  each  consist  of  one  non- zero  term 
as  represented  by  the  left-hand  side  of  Equations  (7)  and 
(8)  .  It  can  be  shown  that  the  terms  of  the  B1  and  # 
matrices  are  of  order  (ka)3  for  small  ka  values.  This  means 
that  we  can  solve  for  the  A-coef ficients  because  we  can 
approximate  the  inverse  of  the  first  large  matrix  of 
Equation  (9) .  By  the  following  development  for  small  ka  we 
can  evaluate  the  coefficients.  We  drop  terms  on  the  order 
of  ka  to  the  sixth  power.  ■ 
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(10) 


I  B'l  1  [i  -B' 

B2  i]  ~[-B2  I 

which  leads  to 

A1 «  C1  -  B‘C2 ;  A2  =  C2  -  B'c1 


Rewriting  this  approximation,  shown  below  in  Equation 
(11) ,  we  then  have  an  explicit  representation  for  the  A- 
coefficients  for  any  mode  of  sphere  one.  We  can  write  a 
similar  equation  for  sphere  two.  Note  that  the  first  term 
is  due  to  its  own  mode  of  vibration;  the  next  term  is  due  to 
sphere  two's  monopole  mode;  the  next  three  terms  are  due  to 
sphere  two's  dipole  modes;  and  the  final  five  terms  are  due 
to  sphere  two's  quadrupole  modes.  The  Kronecker  delta 
symbol  ( S)  equals  one  when  the  n  and  m  values  equal  the 
value  of  the  second  subscript  (i.e.  first  term  exists  when  n 
=  and  m  =  mj  .  In  Equation  (11),  I  have  dropped  the  terms 
on  the  order  of  ( ka )  to  the  sixth  power. 
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E.  SELF  AND  MUTUAL  RADIATION  IMPEDANCE 

The  approximation  for  the  amplitude  coefficients, 
Equation  (11),  plays  an  important  part  of  my  derivation  of 
the  radiation  impedance.  We  will  also  make  use  of  a 
relationship  developed  by  New  and  Eisler  (Ref.  8). 
Radiation  impedance  is  defined  as  the  ratio  of  force  exerted 
by  the  radiator  on  the  median  to  the  velocity  of  the 
radiator.  If  the  velocity  of  the  radiator  is  spatially 
independent,  the  calculation  is  straightforward.  However, 
when  the  velocity  does  have  some  spatial  dependence  then  we 
must  use  the  equation  (*  means  complex  conjugate) 


z„=— 77  iSp(ri)v,(r,)dS,  where  v,  =  V, P(r,) 

\  I /  .  1/  .  surface 

V  J  Y  J  of  the 
sphere 


(12). 


In  applying  Equation  (12),  note  that  the  two  identical 
spheres  may  be  radiating  in  two  distinct  modes.  Using 
Equation  (12)  to  determine  the  radiation  impedance  for 
sphere  one  in  the  presence  of  sphere  two  and  omitting  terms 
on  the  order  of  ( ka )6  we  obtain 


JJIIsEW.afdS, 

rs  /  n,m, 


* 


v: 


1  ‘h*(hnar  (**'*,) 


(13). 
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Reference  9  provides  an  expression  for  the  surface 
integral  of  the  associated  Legendre  polynomials  given  in 
Equation  (13).  Using  this  expression  and  recalling  the 
correct  representation  for  dS  in  spherical  coordinates,  we 
obtain  a  general  form  for  this  surface  integral 

SI  =  2 +|]  (14). 

We  can  evaluate  Equation  (13)  for  several  different 
combinations  of  modes  (e.g.  monopole-monopole,  dipole- 
monopole,  etc...)  and  determine  the  total  radiation  impedance 
for  one  sphere  due  to  both  spheres.  New  and  Eisler 
represent  the  total  radiation  impedance  of  the  first  sphere 
into  the  form  Zrl  =  Z21  +Za*V2/Vt,  where  Zu  is  the  self¬ 

radiation  impedance  and  Z12  is  the  mutual  radiation 
impedance.  For  the  monopole  to  monopole  case  we  obtain  the 
following  expression  for  the  mutual  radiation  impedance. 

t0{“')  <I5) 

Upon  evaluating  several  combinations  of  modes  and 
examining  the  calculated  mutual  radiation  impedance  I  was 
able  to  determine  an  approximate  form  for  the  mutual 
radiation  impedance.  For “  two  identical  spheres  each 
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radiating  in  a  single  mode  of  vibration  the  mutual  radiation 
impedance  is  written  as  (I  omit  terms  on  the  order  of  (ka)5) 


f 


niZ> 'i 

L  a(ni > Px > n2 ’mi>m2)hp (kd) (6n , (f>n ) 

Mn2-n,| 


(16). 


Similarly  the  general  form  for  the  self-radiation 
impedance  for  our  sphere  radiating  in  a  general  mode  can  be 
derived.  Equation  (17)  below  provides  the  general  mode 
self -radiation  impedance  (I  omit  terms  on  the  order  of 
(/ca)5)  . 


-  2m2ip  c  (nt  +  m, ) !  fe.,  ( ka ) 


(17) 


These  two  expressions  provide  a  quick  approximation  for 
the  mutual  and  self  radiation  impedances  between  two  spheres 
each  radiating  in  some  general  mode.  I  will  hereafter  refer 
to  the  mutual  and  self  impedances  calculated  using  Equation 
(16)  and  (17)  as  the  "Modal  Pritchard"  approximation.  Note 
that  the  mutual  radiation  impedance  expression  considers 
orientation  between  the  two  spheres.  The  following  chapter 
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will  investigate  the  utility  of  the  mutual 
impedance  approximation.  Equation  (16). 


radiation 


III.  UTILITY  OF  THE  MODAL  PRITCHARD  APPROXIMATION 

This  chapter  will  investigate  the  utility  of  the  Modal 
Pritchard  Approximation  for  two  types  of  problems .  The 
first  problem  involves  looking  at  the  interaction  between 
two  identical  spheres.  Three  cases  of  two-body  problems  are 
studied:  monopoles,  dipoles,  and  quadrupoles.  The  second 
type  of  problem  involves  one  specific  three-body  where  we 
find  a  limitation  of  our  Modal  Pritchard  approximation  due 
to  its  inability  to  handle  scattering.  Throughout  this 
chapter  I  will  use  a  ka  value  equal  to  one. 

A.  MUTUAL  RADIATION  IMPEDANCE 

1.  Monopole  to  Monopole  Case 

There  are  an  infinite  number  of  combinations  of  mode¬ 
mode  interactions  to  consider,  even  for  the  case  of  two 
identical  spheres.  For  example,  considering  only  up  to  the 
quadrupole  radiating  modes  we  have  9x9=81  combinations  of 
modes  between  two  radiating  spheres  to  examine.  I  will 
limit  the  following  analysis  to  just  three  specific  cases. 

The  first  case  will  be  an  examination  of  the  mutual 
radiation  impedance  between  two  spheres  radiating  as 
monopoles,  also  known  as  the  "breathing  mode".  The  result 
obtained  using  the  "Modal  Pritchard"  approximation  will  be 
compared  to  what  I  will  call  the  "Simple  Pritchard" 
approximation.  Named  after  R.L.  Pritchard,  the  Simple 
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Pritchard  approximation  was  developed  to  provide  an 
expression  for  the  mutual  radiation  impedance  for  two 
cylindrical  sources  placed  on  an  infinite  baffle,  and  was 
further  simplified  by  assuming  a  small  ka  value  (Ref.  10) . 


Z,2=-Real(ZH)i 


(Simple  Pr  itchard  approximation) 


Here  Z12  is  the  mutual  radiation  impedance  between  elements  1 
and  2,  and  ZX1  is  the  self  radiation  impedance  of  element  1. 
Despite  the  ’fact  that  this  approximation  was  made  for 
cylindrical  piston  sources  on  a  theoretical  infinite  baffle, 
other  researchers  have  used  Pritchard's  approximation  for 
study  of  other  shaped  sources  (Ref.  8)  .  In  doing  so,  the 
self  radiation  impedance  is  taken  to  be  the  free-field 
value.  For  the  three-dimensional  case,  Equation  (18)  can  be 
written  using  the  spherical  Hankel  function  as 

L.  =  ReafenhhoM  (19) 


I  will  compare  our  Modal  Pritchard  approximation  with 
the  Simple  Pritchard  approximation  and  note  the  differences. 
In  applying  Equation  (19)  or  Equation  (18) ,  I  calculated  the 
self -radiation  impedance  by  considering  a  sphere  in  free 
field  environment. 

To  gauge  the  quality  of  the  Modal  and  Simple  Pritchard 
approximations,  I  have  chosen  to  compare  the  results  of 
these  methods  against  a  strict  application  of  the  Addition 
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Theorem.  Please  see  Appendix  A  for  the  description  of  my 
application  of  the  Addition  Theorem  along  with  the  computer 
codes  used  to  generate  my  results.  I  choose  to  use  MATLAB 
for  my  calculations  because  of  the  user-friendly  nature  of 
this  engineering  tool .  I  validated  my  Addition  Theorem 
results  with  Fortran  codes  written  by  Professor  Scandrett  of 
the  Naval  Postgraduate  School.  Professor  Scandrett  has 
continually  improved  his  code  over  the  last  eight  years, 
using  it  for  several  research  efforts  (Ref.  11). 

Noting  the  appearance  of  infinity  in  Equations  (4)  and 
(5)  ,  we  are  required  to  truncate  to  some  finite  value  for 
practical  calculations.  Mutual  radiation  impedance 
calculations  using  the  "Addition  Theorem"  method  have 
truncated  this  value  to  six.  Even  with  this  value,  the 
computer  code  must  solve  a  large  (98  by  98)  matrix  of 
equations . 

Figure  (2)  shows  the  calculated  real  part  of  the  mutual 
radiation  impedance  (resistance)  between  two  monopole 
radiating  spheres.  Similarly,  Figure  (3)  shows  the  imaginary 
part  (reactance).  Note  first  of  all  that  the  x-axis  of  each 
plot  is  the  ratio  of  the  center-to-center  distance  between 
the  spheres  to  the  radius  of  the  two  (identical)  spheres. 
Additionally,  the  impedance  values  have  been  scaled  by 
&na2pc. 
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For  both  the  resistance  and  reactance  parts  of  the 
impedance  ,we  observe  that  the  results  of  the  Modal  Pritchard 
approximation  are  in  excellent  agreement  with  results  of  the 
Addition  Theorem  method  throughout  the  plot  ranges.  The 
results  for  the  Simple  Pritchard  approximation  do  not  agree 
with  the  Addition  Theorem  results,  especially  at  smaller 
values  of  (d/a),  which  corresponds  to  closely-spaced 
transducers.  An  analysis  of  the  errors  between  the  Modal 
Pritchard  results  and  the  Addition  Theorem  results,  and 
between  the  Simple  Pritchard  results  and  the  Addition 
Theorem  results,  was  conducted.  Figure  (4)  shows  the 
results  of  this  analysis.  For  ka  =  1,  the  Simple  Pritchard 
approximation  breaks  down  when  the  ratio  between  the  center 
to  center  distance  and  the  radius  is  below  7.5. 


Figure  2.  Monopole  to  Monopole  Resistance 
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Monopole  to  Monopole  Reactance 
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Figure  3 .  Monopole  to  monopole  reactance 


2.  Dipole  to  Dipole  Case 

I  conducted  an  analysis  of  the  dipole-to-dipole  mutual 
radiation  impedance,  similar  to  the  above  analysis  I  took 
the  axes  of  the  dipoles  to  be  oriented  along  the  axis 
joining  the  sphere  centers.  In  terms  of  multipolar 
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components,  this  corresponds  to  the  n=l,  m=0  mode  for  each 
sphere.  For  the  Simple  Pritchard  calculation  I  first 
multiplied  the  real  part  of  the  free- field  self -radiation 
impedance  for  the  hard  sphere  (oscillating  in  a  dipole  mode) 
with  the  appropriate  Hankel  function,  as  represented  in 
Equation  (19) .  This  result  was  then  multiplied  by  three. 
The  reason  for  this  correction  factor  is  that  the 
conventional  Simple  Pritchard  approximation  (Equation  (19)) 
does  not  account  for  any  angular  orientation  of  the  two 
spheres.  The  factor  of  three  correction  was  derived  by 
comparing  the  Simple  Pritchard  Equation  (19)  with  the  full 
Addition  Theorem,  in  a  small  ka  limit,  for  the  dipole 
configuration  examined  (Ref.  12)  .  For  this  particular 
analysis,  I  am  only  considering  one  combination  (n1= 1,  ^=0, 
n2=l ,  and  m,= 0)  .  There  are  a  total  of  nine  (3x3)  different 
combinations  of  dipole  to  dipole  radiation  cases. 

Figures  (5)  and  (6)  show  the  real  and  imaginary  parts 
of  the  mutual  radiation  impedance  for  this  dipole-to-dipole 
case.  Again,  note  the  good  agreement  of  the  Modal  Pritchard 
results  with  the  Addition  Theorem  results.  Notice  how  the 
errors  for  the  Simple  Pritchard  approximation  increase  as 
the  ratio  of  sphere's  spacing  to  the  radius  decreases. 
Figure  (7)  shows  the  errors  between  the  results  for  both  the 
Modal  and  Simple  Pritchard  approximations  as  compared  to  the 
results  for  the  Spherical  Addition  Theorem.  While  the 


errors  for  the  Simple  Pritchard  approximation  are  not  too 
bad  for  large  values  of  d/a,  they  become  significant  when 
d/a  is  less  than  8.0,  for  ka  =  1. 


Dipole  to  Dipole  Resistance 
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Figure  5.  Dipole  to  dipole  resistance.  • 
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Dipoles:  Errors  from  Addition  Theorem 


Figure  (7).  Dipoles,  Error  from  the  Addition  Theorem 


3.  Quadrupole  to  Quadrupole  Case 

As  a  final  example,  I  conducted  an  analysis  of  the 
quadrupole-to-quadrupole  mutual  radiation  impedance.  Two 
linear  quadrupoles  (n=2,  m=0)  were  considered,  with  their 
axes  aligned  with  the  axis  joining  the  sphere  centers.  For 
the  simple  Pritchard  calculation  I  first  multiplied  the  real 
part  of  the  free-field  self-radiation  impedance  for  the  hard 
sphere  with  the  appropriate  Hankel  function,  again  using 
Equation  (19).  In  this  case,  however,  it  can  be  shown  that 
Equation  (19)  needs  to  be  multiplied  by  a  correction  factor 
of  five.  Again,  the  reason  for  this  correction  factor,  as 
in  the  dipole  case,  is  that  the  Simple  Pritchard  formula 
does  not  account  for  any  angular  orientation  of  the  two 
spheres.  For  this  particular  analysis,  I  am  only 
considering  one  combination  (^=2,  m^O,  n2=2,  and  0)  . 
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There  are  a  total  of  twenty- five  different  combinations  of 
quadrupole  to  quadrupole  radiation. 

Figures  (8)  and  (9)  show  the  real  and  imaginary  parts 
of  the  mutual  radiation  impedance  for  this  quadrupole- to- 
quadrupole  case.  We  continue  to  note  the  good  agreement  of 
the  Modal  Pritchard  results  with  .  the  Addition  Theorem 
results.  Notice  again  how  the  error  in  the  Simple  Pritchard 
results  increases  as  the  ratio  between  the  .sphere  spacing 
and  sphere  radius  decreases.  Figure  (10)  shows  the  errors 
between  both  the  Modal  and  Simple  Pritchard  results  as 
compared  to  the  results  using  the  Spherical  Addition 
Theorem.  While  the  errors  for  the  Simple  Pritchard 
approximation  are  not  too  bad  at  large  values  of  d/a,  they 
become  significant  when  d/a  is  less  than  8.0,  for  ka  -  1. 
This  continues  to  show  how  the  Simple  Pritchard 
approximation  breaks  down  as  the  transducers '  spacing 
decreases . 
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Reactance  (scaled)  Resistance  (scaled) 


Quadrupole  to  Quadrupole  Resistance 


Addition  Theorem  -  -  -  a-  -  -  Modal  Pritchard  — * —  Simple  Pritchard 


Figure  (8) .  Quadrupole  to  Quadrupole  Resistance 
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dual  error 


Quadrupoles:  Errors  from  Addition  Theorem 


Figure  (10) .  Quadrupoles,  Errors  from  Addition  Theorem 
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B.  EXAMINATION  WITH  THREE  A  BODY  PROBLEM 

To  further  examine  our  Modal  Pritchard  approximation 
for  the  mutual  radiation  impedance,  we  considered  the 
following  three-body  problem  (see  Figure  11) . 
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In  this  problem  we  will  set  spheres  one  and  two  to 
radiate  as  monopoles.  These  two  transducers  will  maintain 
the  same  physical  dimensions  and  operating  characteristics 
as  before  (i.e.  frequency  =  474  Hz,  sound  speed  of  1490 
m/sec,  radius  =  0.5  m,  ka  =  1,  therefore  the  wavelength 
equals  n)  .  Sphere  three  will  be  considered  acoustically 
hard.  We  examine  the  mutual  radiation  impedance  for  sphere 
one  due  to  the  presence  of  all  three  spheres.  Additionally, 
we  will  examine  the  effect  of  moving  the  third  sphere 
horizontally  away  from  the  axis  of  the  other  radiating 
spheres.  We  will  investigate  the  performance  of  our  Modal 
Pritchard  approximation  for  this  three-body  problem  and 
compare  these  results  with  the  results  using  the  full 
Addition  Theorem.  As  shown  in  Figure  (11),  sphere  three's 
starting  distance  from  the  axis  is  one  quarter  of  a 
wavelength  ( D  -  71/4)  . 

Figure  (12)  shows  a  plot  of  the  mutual  radiation 
resistance,  comparing  the  results  of  our  modal  Pritchard 
approximation  with  those  of  the  Addition  Theorem.  The  x- 
axis  for  this  plot,  as  well  as  Figures  (13)  and  (14),  is  the 
ratio  of  the  distance  ( D )  to  the  radius  of  the  spheres. 
Additionally,  the  results  have  been  scaled  by  Ana  pc. 
Figure  (13)  shows  the  reactance  (or  imaginary)  portion  of 
the  mutual  radiation  impedance.  Figure  (14)  shows  the 
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relative  error  in  the  magnitudes  of  the  results  for  the 
Modal  Pritchard  approximation  against  those  for  the  full 
Addition  Theorem.  Discussions  explaining  my  results  follow 
Figure  (14) . 
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Relative  Error:  Modal  Pritchard  versus  Addition  Theorem 
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The  modal  impedance  approximation  as  I  am  using  it  here 
neglects  scattering.  This  is  why  we  see  the  straight  lines 
for  the  Modal  Pritchard  results  in  Figures  (12)  and  (13)  . 
Because  the  Modal  Pritchard  approximation  does  not  account 
for  the  scattering,  this  provides  an  argument  for  using  the 
full  Addition  Theorem  to  calculate  impedance  when  scattering 
is  important .  Appendix  B  provides  the  method  and  computer 
code  I  used  to  calculate  the  Addition  Theorem  results  for 
this  three-body  problem. 

I  defined  the  relative  error,  as  plotted  in  Figure 
(14),  as  (Modal  Pritchard-Addition  Theorem) /Addition 
Theorem.  Therefore,  dips  in  Figure  (14)  represent 
occurrences  when  the  Addition  Theorem  results  have  a  larger 
magnitude  than  the  magnitude  of  the  Modal  Pritchard 
approximation  results.  Conversely,  peaks  represent 
occurrences  for  which  the  Addition  Theorem  results  are 
smaller  in  magnitude.  Through  the  use  of  the  Addition 
Theorem,  one  can  show  that  these  peaks  and  valleys  result 
from  the  effect  of  scattering. 

The  data  used  to  generate  Figure  (14)  provides  the  D/a 
values  where  these  valleys  and  peaks  occur.  The  first  two 
valleys  occur  at  the  following  D/a  values:  2.57  and  6.07. 
The  first  two  peaks  occur  at  the  following  D/a  values:  4.57 


I  propose  and  will  sketch  out  how  these  valleys  and 
peaks  occur,  due  to  scattering  effects,  by  evaluating  the 
sum  of  the  pressure  waves  arriving  at  sphere  one,  I  will 
limit  this  evaluation  by  only  looking  at  the  monopole  case. 
Since  the  Addition  Theorem  correctly  handles  scattering,  the 
Addition  Theorem  results  provide  local  maximum  values  when 
the  scattering  waves  from  sphere  three  are  in  phase  with  the 
direct  path  waves  from  sphere  two.  Similarly,  the  Addition 
Theorem  results  have  local  minimum  values  when  the 
scattering  waves  are  180°  out  of  phase  with  the  direct 
path's  phase.  Our  Modal  Pritchard  approximation  does  not 
handle  the  scattering  effects,  note  the  flat  plots  in 
Figures  (12)  and  (13)  for  the  Modal  Pritchard  results. 

Using  the  Addition  Theorem  and  considering  only  the 
monopole  case,  the  pressure  waves  are  of  the  form  of 
outgoing  spherical  Hankel  function  as  shown  below. 

Ib'{kr)=  j0(kr)fi2’(kd:i)  (20) 

An  examination  of  the  direct  path  wave  from  sphere  two 
to  sphere  one  and  of  the  reflected  wave  from  sphere  two, 
reflected  off  sphere  three  towards  sphere  one,  was  completed 
by  Professor  Scandrett  (Ref.  12)  .  The  expression  for  the 
sum  of  these  waves  at  sphere  one  is  provided  below.  The 
first  term  in  the  bracket  is  due  to  the  direct  path  and  the 
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second  term  is  due  to  the  reflected  (scattered)  wave.  The  A 
is  some  arbitrary  amplitude  coefficient. 


Aj„(k,a) 


A  plot  of  Equation  (21)  is  provided  in  Figure  (15)  in 
terms  of  D/a  versus  magnitude.  The  peaks  of  Figure  (15) 
match  the  valleys  of  Figure  (14)..  Similarly,  local  minimums 
of  Figure  (15)  match  the  peaks  of  Figure  (14).  The  first 
two  peaks  of  Figure  (15)  are  at  D/a  values  of  2.57  and  6.07 
and  the  first  two  local  minimums  are  at  D/a  values  of  4.57 
and  7.82.  The  maximum  values  occur  when  the  scattering 
waves  are  in  phase  with  the  direct  path  waves.  The  minimum 
values  of  Equation  (21)  occur  when  the  scattering  wave  is 
180  degrees  out  of  phase  with  the  direct  path. 


Figure  (15) :  Plot  of  Equation  (21) 
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Figure  (16)  provides  an  additional  drawing  of  this 
three-body  problem  and  shows  the  distance  parameters  (d12  and 
d13)  used  in  Equation  (21). 


Again,  the  full  Spherical  Addition  Theorem  is 
recommended  over  our  Modal  Pritchard  approximation  when 
scattering  is  an  important  consideration  as  in  this  three- 
body  problem. 


37 


(THIS  PAGE  INTENTIONALLY  LEFT  BLANK) 


IV.  NUMERICAL  MODELING  WITH  ATILA 

This  chapter  describes  the  ATILA  finite-element  model 
used  to  calculate  the  scattering  of  acoustic  waves  from  an 
elastic  sphere  in  water.  The  goal  is  to  compute  the  so- 
called  T-matrix,  which  relates  the  scattered  waves  to .  the 
incident  waves,  in  an  orthogonal  (spherical  harmonic)  basis 
set  (Ref.  11). 

A.  MODEL  INTRODUCTION 

A  few  comments  about  ATILA  are  appropriate  at  this 
time.  As  mentioned  in  the  introduction,  ATILA  is  a  French 
designed  finite  element  code  specifically  developed  to  model 
sonar  transducers.  ATILA,  written  mostly  in  standard 
FORTRAN  77,  is  the  result  of  many  years  of  research  done  at 
the  Insitut  Superieur  d'Electronigue  du  Nord.  It  is  quite 
large,  consisting  of  over  200,000  lines  of  code  (Ref.  13). 

This  tool  has  been  used  by  NPS  faculty  and  masters 
thesis  students  for  several  years  and  has  been  updated  at 
least  three  times.  Recent  studies  have  helped  to  identify 
needed  modifications  to  the  code,  such  as  increasing  the 
maximum  allowable  number  of  degrees  of  freedom,  and 
modification  of  the  numerical  techniques  used  to  handle  the 
radiation  boundaries.  For  example,  during  the  calculation 
of  the  entries  for  the  transition  matrix  (T-matrix)  there 
were  noted  errors  when  compared  to  answers  derived  from  a 
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known  analytical  solution  to  the  problem  (Ref.  11).  The 
ATILA  version  I  used  (5.1.1)  did  increase  the  maximum 
allowable  degrees  of  freedom,  but  has  yet  to  address  the 
radiation  boundary  problem. 

B.  ATILA  AND  MESH  DESCRIPTION 
' 1.  About  the  Model 

A  distinguishing  characteristic  of  ATILA  is  the  use  of 
isoparametric  elements  for  both  the  shell  elements  and  the 
fluid  elements.  Isoparametric  elements  are  elements  that 
use  the  same  shape  function  for  geometric  mapping  as  well  as 
field,  calculations  (Ref.  4).  The  major  advantage  of 
isoparametric  elements  is  that  they  lend  themselves  well  to 
numerical  integration.  This  is  because  of  the  nature  of 
isoparametric  elements.  These  elements  are  defined  in  the 
natural  domain  to  be  normalized,  thus  making  numerical 
integration  much  easier  to  apply  (Ref.  4).  Also,  when 
isoparametric  elements  are  used,  fewer  nodes  are  needed  to 
define  the  mesh  than  without  such  elements. 

Other  numerical  techniques  for  solving  the  acoustic 
equations  have  been  used  by  U.S.  Navy  labs.  One  method  used 
for  many  years  at  the  Naval  Ocean  Systems  Center  is  called 
CHIEF  (Combined  Helmholtz  Integral  Equation  Formulation) 
(Ref.  2).  CHIEF,  which  models  acoustic  radiation  from 
bodies  with  arbitrary  shapes  does  not  use  isoparametric 
elements.  Thus  we  should  expect  a  higher  degree  of  accuracy 
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with  ATILA  over  CHIEF.  One  advantage  to  CHIEF  is  that  the 
radiation  boundary  condition  is  "built-in"  to  the  boundary 
element  by  its  use  of  the  Helmholtz  Integral  Equation. 
Currently,  this  provides  an  advantage  over  ATILA,  which 
employs  so-called  "monopolar"  and  "dipolar"  "dampers". 
These  dampers  can  not  represent  a  perfectly  absorbing 
radiation  boundary  for  a  higher  order  multipolar  field. 

The  basic  set  of  equations  that  allow  for  a  large 
number  of  analysis  to  be  performed  include  elasticity  in  the 
structure,  the  Helmholtz  equation  in  the  fluid,  and 
Poisson's  equation  in  the  elastic  or  piezoelectric  material. 
Primary  variables  found  using  ATILA  include  the  displacement 

field  U  in  the  whole  structure,  the  electric  potential  &  in 

the  piezoelectric  material,  and  the  pressure  P  in  the  fluid. 
In  matrix  form  the  equations  are  as  follows  (Ref.  13) . 


(LsrJ -o)2\m  ]) 
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(22) 


Vectors  in  the  above  equations  are: 

U=  nodal  values  of  the  components  of  the  displacement  field, 
0=  nodal  values  of  the  electric  potential, 
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P=  nodal  values  of  the  pressure  field, 

F=  nodal  values  of  the  applied  forces, 

g=  nodal  values  of  the  electric  charges,  and 

nodal  values  of  the  integrated  normal  derivative  of  the  pressure  on  the 
surface  boundary  S. 

Matrices  in  Equation  (22)  include 

Kuu  =  stiffness  matrix, 

Ku<p  =  piezoelectric  matrix, 

K^,  =  dielectric  matrix, 

M  =  consistent  mass  matrix, 

H  =  fluid  (pseudo-)  stiffness  matrix, 

Mi  =  consistent  (pseudo-)  fluid  mass  matrix, 

L  =  coupling  matrix  at  the  fluid  structure  interface,  and 

O  =  zero  matrix. 

In  addition, 

ct)=  Angular  frequency, 
p = Fluid  density, 
c  =  fluid  sound  speed,  and 
T  =  means  transposition  of  a  matrix. 

2 .  Three-dimensional  Fluid  Mesh 

To  run  this  and  other  finite  element  models  you  need  to 
split  the  region  under  study  into  elements;  the  ATILA  code 
has  an  extensive  library  of  elements  available.  These 
elements  are  distinguished  by  a  set  ordering  of  nodes  and 
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the  node  coordinates.  Figures  16  provides  a  Mercator 
projection  of  the  eight  "super-elements"  used  for  one  of  the 
several  layers  under  study.  These  super  elements  are  the 
start  of  the  mesh  generation.  The  solid  lines  and  nodes 
(dots  with  numbers)  define  the  super-elements.  The  dotted 
lines  show  how  the  super- elements  are  sub-divided.  Laid  out 
3-dimensnionally,  Figure  17  represents  the  super-elements 
that  make  up  the  spherical  transducer  under  study,  a 
spherical  shell  with  a  radius  of  0.5  m.  Using  a  finite 
element  mesh  generator  called  MOSAIQUE,  I  built  a  complex 
mesh  of  the  entire  field  using  the  super-elements,  using 
special  instructions  for  subdividing  the  elements  into  a 
large  number  of  elements  and  nodes  that  will  define  each 
layer  of  the  field.  Figure  18  provides  a  detailed  look  of 
the  field  under  study  (the  fluid  outer  boundary  is  at  5 
meters) . 

There  are  several  layers  that  make  up  the  entire  mesh 
some  of  which  are  full-layers  and  some  are  mid- layers.  Each 
full-layer  consists  of  194  nodes  of  which  18  are  nodes  of 
the  super-elements.  Mid-layers  have  a  total  of  62  nodes. 
Starting  from  the  innermost  radius  and  working  outward,  the 
complete  three-dimensional  mesh  consists  of  the  following 
layers.  First  is  the  shell  layer  that  consists  of  194 
nodes.  This  is  the  actual  spherical  shell  transducer  that 
is  being  analyzed.  Next,  another  194  nodes  are  used  to 
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describe  the  interface  elements  used  to  model  the  shell- to- 
water  interface.  This  layer  shares  the  same  coordinate 
locations  as  the  shell's  nodes.  Following  this  are  ten 
fluid  layers,  the  first  two  of  which  are  0.25  meter  thick 
and  the  remaining  are  1.0  meter  thick.  Each  fluid  layer 
consists  of  a  full-layer  and  a  mid-layer.  Finally 
completing  this  three-dimensional  mesh  we  have  the  radiation 
boundary  elements.  These  last  elements  are  used  to 
prescribe  monopolar-only  or  monopolar  and  dipolar  radiation 
boundary  conditions  and  are  used  to  terminate  the  mesh. 


18  18  18  18 

Figure  17.  The  Super-Element  (Macerator  projection) 
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Figure  19.  The  complete  Mesh 
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V.  NUMBERICAL  MODELING  WITH  ATILA:  RESULTS  AND  DISCUSSION 

A.  OBJECTIVES  AND  THE  T-MATRIX 

My  objectives  for  the  numerical  modeling  portion  of 
this  thesis  were  to  investigate  the  finite  element  mesh 
mentioned  earlier,  simplify  the  process  of  calculating  the 
Transition  Matrix  (or  T-Matrix) ,  and  provide  a  new  method  of 
evaluating  the  off-diagonal  elements  of  the  T-matrix.  The 
finite  element  code  is  used  to  numerically  calculate  the 
scattered  pressures  from  an  object  tinder  study.  We  then  use 
these  results  to  calculate  the  object's  specific  T-Matrix. 

In  this  chapter,  I  will  document  the  results  of  the  T- 
matrix  calculation  and  compare  with  previous  numerical 
modeling  efforts.  Several  references  mention  and  thoroughly 
discuss  the  definition,  usage,  and  calculation  of  the  T- 
Matrix  (Ref.  3  and  5)  .  I  will  simply  mention  here  in  the 
text  that  the  T-Matrix  describes  the  scattering  properties 
of  a  specific  body.  It  does  this  by  using  a  discrete  basis 
set  spherical  harmonics  in  our  case,  and  is  used  to  show  how 
the  incident  pressure  waves  translate  to  scattered  waves 
from  the  body.  Please  see  Appendix  C  for  more  information 
about  the  T-Matrix  and  how  I  calculated  the  elements  of  the 
T-Matrix  for  our  radiating  body. 
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B.  RESULTS 

In  order  to  gauge  the  accuracy  of  our  T-Matrix 
calculation  we  compare  our  answers  with  an .  analytically 
derived  equation.  An  expanded  discussion  of  the  derivation 
of  the  analytical  solution  can  be  found  in  a  Master's  Thesis 
by  Ruiz  (Ref.  5  and  references  therein).  For  a  thin 
spherical  shell  the  T-Matrix  is  diagonal,  and  are  given  by 


iz,hjka)  +  p,c,h,(ka) 


(23) 


where 
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,/s,  (0-2nf ,  An-n{n  + 1) 


cr  '  \i  A  (!">') 

v  =  Poisson's  ratio  (0.33);  /z  =  shell  thickness  (1  cm) 
cp  =  shell's  plate  velocity 
E  =  Young’s  modulus  (0.125 xlO12) 
a  =  spherical  shell  radius  (0.5  meter) 

p  =  shell's  density  (7500  kg  fm3);f  =  frequency  (474  #z) 
p  and  cf  =  fluid  density  (1000  kg  fm3)  and  speed  (1490  m  /  5). 


The  calculation  for  this  analytical  solution  was 
performed  using  MATLAB  and  the  computer  program  can  be  found 
at  the  end  of  Appendix  C.  In  the  following  discussion,  when 
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referring  to  errors,  I  am  comparing  the  calculated  values 
with  the  analytical  solutions  from  Equation  (23) . 

Using  the  thin  spherical  shell  parameters  listed  above 
and  the  MATLAB  computer  script  provided  in  Appendix  C,  I 
calculated  the  analytical  T-Matrix  diagonal  elements  (the 
off-diagonal  elements  are  zero  for  a  spherically  symmetric 
scatter) .  These  diagonal  elements  are  provided  in  Table  I 
below.  A  value  of  ka  =  1  was  used  (frequency  =  474  Hz) . 


ANALYTICAL  DIAGONAL  T-MATRIX  ENTRIES 

ELEMENT 

REAL 

IMAGINARY 

MAGNITUDE 

PHASE 

PART 

PART 

(DEGREES) 

Tn  =  Roo 

-1.3465E-02 

1.1525E-01 

1.1604E-01 

96.6635 

T22  =  Ri-i 

-6.0255E-03 

7.7390E-02 

7.7624E-02 

94.4520 

T33  =  Rio 

-6.0255E-03 

7.7390E-02 

7.7624E-02 

94.4520 

T44  =  Rn 

-6.0255E-03 

7.7390E-02 

7.7624E-02 

94.4520 

T  55  -  R2-2 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

T  66  =  R2-I 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

T77  =  R20 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

T88  "  ^21 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

Tgg  =  R22 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

Table  I.  Analytical  Diagonal  T-Matrix  Entries. 


I  then  used  the  finite  element  code  (ATILA)  to  model 
the  scattering  pressure  field  generated  by  the  thin 
spherical  shell.  With  the  scattered  pressure  results,  I 
calculated  the  T-Matrix  entries  column  by  column.  These 
results  are  provided  below  in  Table  II . 
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Numerical  Diagonal  T-Matrix  Elements 

Element 

Real 

Imaginary 

Magnitude 

Phase 

Part 

Part 

(Degrees) 

T„ 

-1.6107E-02 

1.2894E-01 

1.2994E-01 

97.1202 

T22 

-8.0999E-03 

7.9943E-02 

8.0353E-02 

95.7855 

T33 

-7.7968E-03 

7.9777E-02 

8.0157E-02 

95.5819 

T44 

-8.080 IE-03 

7.9946E-02 

8.0353E-02 

95.7713 

t55 

-2.1393E-03 

-1.5865E-02 

1.6009E-02 

-97.6794 

T66 

-2.2524E-03 

-1.5923E-02 

1.6082E-02 

-98.0513 

T77 

-2.1536E-03 

-1.5930E-02 

1.6075E-02 

-97.6994 

T88 

-2.250 IE-03 

-1.5925E-02 

1.6084E-02 

-98.0420 

T99 

-2.1498E-03 

-1.5860E-02 

1.6005E-02 

-97.7193 

Table  II.  Numerical  Diagonal  T-Matrix  Elements. 


Table  III,  next  page,  provides  an  error  analysis  for 
the  current  results.  Unfortunately,  the  monopole  result 
(i.e.  the  first  diagonal  element)  is  much  worse  than 
previous  results.  All  other  diagonal  elements  show  some 
improvement  over  previous  work,  between  6  and  56  percent 
improvement  in  their  magnitudes.  But  the  most  important 
indicator  is  the  monopole  element;  therefore,  our  new  mesh 
does  not  provide  a  significant  improvement  to  the  T-Matrix 
calculation.  As  the  ATILA  code  is  further  improved  we  hope 
that  these  errors  will  be  corrected. 
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Error  Comparison  Current  vs.  Previous 

Element 

Magnitude 

Phase 

Previous 

Previous 

Magnitude 

Phase 

Error 

Error 

Mag.  Error 

Phase  Error 

Improve 

(percent) 

(Degrees) 

(Degrees) 

(percent) 

(Degrees) 

Tn 

11.984 

0.457 

-0.360 

-1216.94 

-0.10 

T22 

3.515 

1.333 

6.490 

0.329 

45.84 

-1.00 

T33 

3.263 

1.130 

7.330 

-0.500 

55.48 

-0.63 

T44 

3.516 

1.319 

6.470 

0.314 

45.66 

-1.01 

T55 

-35.588 

-6.255 

-37.800 

-9.186 

5.85 

2.93 

Teg 

-35.296 

-6.627 

-38.000 

-9.006 

7.12 

2.38 

t77 

-35.324 

-6.275 

-37.700 

-7.823 

6.30 

1.55 

t88 

-35.288 

-6.618 

-8.986 

7.14 

2.37 

T99 

-35.605 

-6.295 

-37.900 

-9.106 

6.05 

2.81 

Table  III.  Error  Comparison  Current  vs.  Previous. 

The  T-Matrix  must  be  diagonal  due  to  the  symmetry  of 
the  object.  That  means  that  the  off-diagonal  elements  must 
be  zero.  Previous  work  documented  eight  significant  off- 
diagonal  elements.  Current  results  reveal  these  same  eight 
elements  are  the  most  significant  off-diagonals.  However, 
the  new  results  show  improvement,  meaning  the  new  off- 
diagonal  values  are  much  closer  to  zero  than  previous 
calculations.  Along  with  matching  the  analytical  diagonal 
elements,  making  all  off-diagonal  elements  closer  to  zero  is 
another  indication  of  how  well  our  finite  element  method  is 
working.  Table  IV  offers  a  comparison  of  this  improvement. 
The  last  column  of  this  table  shows  the  percentage  closer  to 
zero  the  new  off-diagonal  magnitudes  are  compared  to  the 
previous  work. 
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SIGNIFICANT  OFF-DIAGONALS 


Element 

Real 

Part 

Imaginary 

Part 

Improvement 
(percent  closer  to  zero) 

t7, 

3.21650E-04 

6.87400E-04 

18.96 

t62 

-9.89950E-08 

-1.78470E-07 

99.12 

Tg4 

-3.02040E-18 

8.25020E-19 

100.00 

T95 

1.18910E-07 

-9.80070E-07 

95.46 

T26 

3.69920E-06 

1.98190E-06 

62.58 

Tn 

-5.87890E-04 

1.30000E-03 

79.74 

T48 

9.56520E-07 

2.56450E-07 

98.07 

T59 

2.68580E-04 

1.18460E-04 

97.69 

Table  IV.  Significant  Off-Diagonals. 


C.  FURTHER  ANALYSIS  AND  DISCUSSION 

Some  Matrix  algebra  is  employed  to  further  characterize 
how  well  the  T-Matrix  off-diagonals  were  calculated.  The 
following  equation  provides  a  normalization  method  used  to 
characterize  the  off-diagonal  elements  relation  to  one. 
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By  multiplying  both  sides  of  the  T-Matrix  with  the 
diagonal  matrix  as  described  in  Equation  (24)  I  will  obtain 
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a  matrix  with  one's  along  the  diagonal  and  relative  errors 
in  the  off-diagonal  elements.  By  examining  the  magnitudes 
of  the  most  significant  off-diagonal  elements,  I  can  quickly 
gauge  how  "well"  I  calculated  the  T-Matrix.  Again  we  see 
the  same  eight  off-diagonals  show  the  most  significant 
levels  with  respect  to  the  normalization.  Although  they  are 
the  most  significant  they  are  small.  This  method  of 
analyzing  the  T-Matrix' s  off-diagonals  provides  a  "quick- 
check"  as  to  how  well  our  numerical  analysis  is  working. 
Table  V  provides  our  results  for  this  analysis.  Notice  how 
small  the  magnitudes  of  these  elements.  All  other  elements 
have  magnitudes  on  the  order  of  element  TM  and  smaller. 


T-MATRIX  NORM/ 

(most  significa: 

lLEZATTON  ANALYSIS 

NT  OFF-DIAGONALS) 

Element 

Real 

Imaginary 

Magnitude 

Part 

Part 

t7, 

6.96163E-03 

1.50759E-02 

1.6606E-02 

t62 

-2.655 19E-06 

-5.01816E-06 

5.6773E-06 

T84 

-8.44548E-17 

2.12802E-17 

8.7095E-17 

T95 

5.96808E-05 

1.55646E-05 

6.1677E-05 

t26 

1.01798E-04 

5.71578E-05 

1.1675E-04 

Tit 

-1.30097E-02 

2.89475E-02 

3.1737E-02 

t48 

2.64606E-05 

7.65933E-06 

2.7547E-05 

T59 

-9.58209E-03 

1.56365E-02 

1.8339E-02 

Table  V.  T-Matrix  Normalization  Analysis. 

! 
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VI.  CONCLUSION  AND  FOLLOW-ON  PROPOSALS 

A.  CONCLUSION 

i 

I  have  derived  an  approximation  for  the  self  and  mutual 
radiation  impedances  for  spherical  transducers  of  small  ka 
values  by  considering  two  spheres  in  isolation.  In  this 
thesis,  I  have  decided  to  call  this  approximation  "Modal 
Pritchard" .  This  approximation  was  tested  for  a  general 
sphere  of  radius  a= 0.5  meters,  frequency  of  474  Hz,  in  water 
with  sound  speed  of  1490  m/sec.  Two  general  problems  were 
tested:  (1)  two  bodies  radiating  in  three  modes  (monopoles, 
aligned  dipoles,  and  aligned- linear  quadrupoles)  and  (2)  a 
three-body  problem  in  which  two  bodies  radiated  in  an  array 
and  the  third  acts  as  a  moving  scatterer.  In  all  the  two- 
body  problems,  the  Modal  Pritchard  approximation  provides 
mutual  radiation  impedance  results  close,  if  not  exactly, 
matching  the  full  Addition  Theorem  results.  The  same  cannot 
be  said  for  the  three-body  problem  where  scattering  is 
important.  Therefore,  when  scattering  is  important  we 
recommend  using  the  full  Addition  Theorem.  In  some  cases, 
our  Modal  Pritchard  approximation  may  be  a  useful  tool  in 
the  design  phase  of  low-frequency  active  sonars. 

The  numerical  modeling  portions  of  this  thesis  have 
updated  and  streamlined  the  process  of  calculating  the 
scattering  characteristics  and  T-matrix  for  an  object  under 
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study.  Using  a  more  refined  mesh,  only  a  partial 
improvement  has  been  realized  for  the  T-matrix  calculation 
for  a  thin  spherical- shell  radiator.  Unfortunately,  current 
work  showed  worse  errors  for  the  monopole  element  of  the  T- 
Matrix . 

We  also  noted  some  improvements  in  the  calculation  of 
the  off-diagonal  elements.  Most  off-diagonal  elements 
showed  an  improvement  of  between  63  to  100  percent  closer  to 
zero  as  compared  to  previous  work.  I  also  provided  a  quick 
and  easy  means  for  determining  how  well  we  calculated  the  T- 
matrix  off-diagonal  elements  by  "normalizing"  the  matrix. 
With  this  process,  we  can  quickly  gauge  the  off-diagonal 
elements.  AS  a  result,  seven  off-diagonal  elements  were 
between  0.02  and  0.000006  as  compared  to  the  number  one. 
All  remaining  off-diagonal  elements  (65)  were  smaller  than 
9.0*10'17  (most  significantly  smaller).  Current  work  shows 
our  T-matrix  calculation  was  worse  than  previous  work  for 
the  monopole  element,  but  showed  some  improvements  for  other 
elements  of  the  T-matrix.  There  is  still  room  for  further 
improvement  which  may  be  realized  when  ATILA  radiation 
boundary  problems  are  addressed. 


B.  FOLLOW-ON  PROPOSALS 


Further  work  with  the  "Modal  Approximation"  to  the 
mutual  radiation  impedance  include  studying  other  problems 
such  as  arrays  of  multiple  transducers  and  calculating  the 
pressure  field  generated  from  their  radiation.  This 
investigation  does  not  need  to  be  limited  to  comparisons 
with  computer  generated  results  from  the  Addition  Theorem. 
Follow-on  work  could  use  this  Modal  Pritchard  approximation 
to  compare  with  laboratory  results  of  real  transducers . 
With  meticulous  data  gathering,  one  could  study  how  well 
this  approximation  does  in  the  "real"  world. 

As  mentioned  earlier,  there  continues  to  reside 
problems  with  the  radiation  boundary  treatment  in  ATILA. 
This  problem  was  to  have  been  addressed  with  an  even  newer 
version  of  ATILA  to  include  a  tool  called  EQI.  In  reference 
to  the  second  portion  of  this  thesis,  it  is  hoped  that  this 
modification  will  further  improve  the  T-matrix  calculation. 
Currently,  about  half  of  our  T-matrix  diagonal  elements  show 
magnitude  errors  below  7.3  percent  and  the  rest  are  about  38 
percent  in  error.  Current  phase  errors  are  below  9  percent. 
Further  analysis  using  the  modified  ATILA  with  EQI  may  reach 
our  goal  of  relative  errors  below  one  percent.  Future  work 
could  also  use  other  refined  meshes  that  may  also  provide 
further  improvements. 
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APPENDIX  A:  ADDITION  THEOREM  CALCULATION  OF  MUTUAL  RADIATION 

IMPEDANCE 

Here  I  will  explain  how  I  used  the  addition  theorem  to 
calculate  the  mutual  radiation  impedance  for  two  identically 
hard  spherical  shells.  I  will  refer  to  some  equations  in 
the  text  of  this  thesis  throughout  this  explanation  and  will 
provide  my  MATLAB  program  along  with  some  user  defined 
functions. 

Recall  the  application  of  the  Addition  Theorem  to  write 
the  pressure  from  one  sphere  relative  to  the  coordinate 
system  of  the  other  sphere.  Equations  (4)  and  (5)  of  the 
text.  The  application  of  the  boundary  conditions  as  written 
in  Equation  (6)  brings  us  to  Equations  (8)  and  (9).  I  will 
rewrite  the  last  two  equations  here. 

°°  A  J  t2+rii  ’  t  \ 

X  X  AL;?‘  X  a  vh'Pi'ti'nU’Sir 
0  12hnmPlf-n,\ 

pAst~m\ 

hp  ikd 

00  A  J  Vn*  ’  /  \ 

X  X  A  l5i-r2— -■  X  a  vh’Pi’ti'nh'Sir 

trO  sr-ti  tlS'hn(^)p^trn2\ 

Pt^srmJi 

hpikd  )Q?fm'-in-dn,n+<t>l) 
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These  two  equations  are  put  into  matrix  form  as  given 
by  Equation  (10)  .  Notice  the  infinity  sign  in  the 
expression  above.  For  practical  purposes  we  have  to 
truncate  this  infinity  to  some  finite  number  K.  Using  K  = 
6,  the  left  hand  side  of  Equation  (10)  becomes  a  large  98  by 
98  matrix  multiplying  a  98  element  column  vector  for  the 
left  hand  side  and  the  right  hand  side  is  98 -element  column 
vector  with  only  two  non-zero  entries.  We  solve  for  the 
amplitudes  in  MATLAB  by  using  a  forward  divide  as  follows. 


With  K  =  6  the  matrices  have  the  following  sizes, 


Now  we  have  solved  the  unknown  A-coef f icients  and  are 
ready  to  use  these  results  to  find  the  mutual  radiation 
impedance . 

Using  Equation  (13)  with  our  solved  A-coef ficient 
values,  we  can  easily  solve  for  the  radiation  impedance. 
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given  that  the  total  radiation  impedance  can  be  written  as 
Zrl  =  Z21  +Z12*V2/V1.  I  first  solve  for  the  total  radiation 
impedance  by  setting  both  V' s  equal  to  one.  Then  by  setting 
V2  =  0  I  obtain  the  self-radiation  impedance  (Z21)  result.  I 
then  subtract  Z12  from  Zrl  to  obtain  the  mutual -radiation 
impedance  ( Z12 )  result. 

Here  are  the  MATLAB  computer  scripts  used  to  generate 
these  results  along  with  all  supporting  user  defined 
functions  written  for  this  analysis. 


MAIN  MATLAB  PROGRAM 
%  Joseph  L.  Day 

%  Mutual  impedance  using  the  Addition  Theorem 
%  2  spheres  along  z-axis,  with  top  sphere  moving  away 
%  Last  updated  July  27,  1999 
format  long 

K=6;  %  Truncated  integer  vice  infinity 

K2=(K+1)A2;  %  Useful  number,  it  is  the  size  of  many  vectors 
k=2;  %  Acoustic  wave  number. 

radius=0 . 5;ka=k*radius;  %  Both  spheres  with  same  radius. 
theta=0;phi=pi;  %  Angle  from  2's  origin  to  l's 

nl=0;  ml=0;  n2=0;  m2=0;  %  Defines  the  mode  for  each  of  the 

two  spheres 

positional;  %  Positionl  and  Position2  are  integers  change 
position2=l;  %  as  follows:  nl=0,ml=0  then  positionl=l; 

%  nl=l,ml=-l  then  positionl=2;  and 
%  nl=l,ml=0  then  positionl=3, . .and  so  on. 
V1=1;V2'=0;  %  Spheres  velocity  amplitudes 

%  Both  equal  to  one  gives  Z (total) 

%  With  Vl=l  and  V2=0  gives  Zll 
speed=1490;  %  speed  of  sound  in  the  water. 

rho=1000;  %  density  of  the  water. 

templ=eye (K2 ) ; temp2=zeros (K2 , K2 ) ; 

Big= [tempi  temp2;temp2  tempi];  %  Matrix  allocation 
count=l;  %  Initialize  a  counting  process 

for  d=2.5:.25:24  %  Distance  with  0.25  meter  increment 
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countl=l ; count2=l ;  %  Counting  process. 

Btempld=zeros ( 1 ,  K2 ) ; 

Btemp2d=zeros (1, K2 ) ; 

for  t=0 : 6  . 

for  s=-t: t 

ptempla=abs ( t-nl) ;ptemplb=abs (s-ml) ; 
ptemplc=max(ptempla,ptemplb) ; 
ptemp2a=abs ( t-n2 ) ; ptemp2b=abs ( s-m2 ) ; 
ptemp2c=max(ptemp2a,ptemp2b) ; 
count3=l ; count4=l ; 
for  pl-ptemplc : 2 : t+nl 
M=s-ml; 
if  pl<0 

P=  (-1) Aabs (M) *fact (pl-abs (M) )/ 
fact (pl+abs (M) ) *legendre (pi, cos (theta) ) ; 

else 

P=legendre (pi, cos ( theta) ) ; 

end 

Btempla=a(nl,pl,t,ml,s)*shank2(pl,d); 

Btemplb=P(l+abs(M))*exp(i*M*phi); 

Btemplc (1, count3) =Btempla*Btemplb; 
count3=count3+l; 

end 

Btempld(l,  count l)  =sum( Btemplc)  *jprime  (nl, ka)  /hprime  (nl, ka)  ; 
countl=countl+l; 
for  p2=ptemp2c : 2 : t+n2 
M=s-m2 ; 
if  p2<0 

P= (-1) Aabs (M) *fact (p2-abs (M) ) / 
fact (p2+abs (M) ) *legendre(p2, cos (theta-pi) ) ; 

else 

P=legendre (p2 , cos ( theta-pi )) ; 

end 

Btemp2a=a(n2,p2,  t,m2,s)  *shank2  (p2,d)  ; 
Btemp2b=P(l+abs(M))*exp(i*M*(phi+pi)); 

Btemp2c (1, count4) =Btemp2a*Btemp2b; 
count4=count4+l ; 

end 

Btemp2d ( 1 , count 2 ) =sum(Btemp2c) * jprime (n2 , ka) /hprime (n2 , ka) ; 
count2=count2+l ; 

end 

end 
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Big(positionl,K2+l:2*K2)=Btempld;  %Placing  values 
Big (K2+position2 , 1 :K2) =Btemp2d;  %in  the  proper  row. 
C (2*K2 , 1) =zeros; 

C(positionl,l)=-i*rho*speed*Vl/hprime(nl,ka); 

C (position2+K2 , 1) =-i*rho*speed*V2/hprime (n2 , ka) ; 
A=Big\C; 


Temporary=sum(conj  (A(K2+1 : 2*K2 , 1) )  '  .  *Btempld*sbessj  (nl,ka) )  * 
hprime(nl,ka) /jprime(nl,ka)  ; 

Zrl=4*pi*radiusA2/Vl* (A(l, 1) *shank2 (nl, ka) + Temporary) ; 
Zrvec tor (count) =Zrl ; 
count=count+l ; 

end 

Zrlscaled=Zrvector/ (4*pi*radiusA2*rho*speed) ; 

%  This  last  value  is  the  final  answer  and  has  been 
%  scaled  for  plotting  purposes. 

%  Also  for  plotting  purposes  the  final  answer  is  calculated 
%  for  ranges  (d/a)  from  2.5  to  24  in  increments  of  0.25. 

%  All  plots  were  done  using  Microsoft  Excel,  since  that 
%  program  works  better  with  Microsoft  Word  for  cutting  and 
%  pasting  graphs. 
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Supporting  user  defined  functions : 
function  A=a (s, t, r, u,m) 

%  computes  the  function  alS/tjr/U/m)  of  the  Addition  Theorem 
%  as  written  in  King  and  Van  Buren 

%  Written  by  Joseph  L.  Day 
numl= (2*s+l) * (2*t+l) *fact (s-u) *fact ( t- 
m+u) *fact (r+m) *fact ( (s+t+r) /2 ) ; 
deml=fact( (r+t-s) /2) *fact ( (r+s-t) /2) * 
fact ( (s+t-r) /2) *fact (s+t+r+1) ; 
terml=numl  /deml ;  * 

wmin=0 . 5*max ( [r-s-t, s-r-t-2*m+2*u, t-s-r+2*u] ) ; 
wmax=0 . 5*min ( [s+t-r, r+t-s-2*m+2*u, r+s-t+2*u] ) ; 
sum=0; 

for  w=wmin:wmax 

sum=sum+ (-1) Aw*bc (s+t-r, (s+t-r) /2+w) * 

bc(t+r-s, (t+r-s) /2+m-u+w) *bc (s+r-t, (s+r-t) /2-u+w) ; 

end 

A=real (iA (s+t-r) ) *terml*sum;  %  Output  from  this  function 


function  y=bc(n,m) 

%  Computes  the  binomial  coefficient 

%  Written  by  Joseph  L.  Day 
if  m<0 
y=0; 

elseif  n-m<0 
y=0; 
else 

y=prod (1 :n) / (prod (1 :m) *prod (1 :n-m) ) ; 

end 


fiinction  y=fact(n) 

%  Computes  the  factorial  of  n 
%  Written  by  Joseph  L.  Day 
y=prod(l:n); 
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function  hn=shank2 (n,x) 

%  Computes  the  Spherical  Hankel  function  of  the  second  kind 

%  Written  by  Joseph  L.  Day 

hn=sqrt (pi/ (2*x) ) *besselh(n+ . 5 , 2 ,x) ; 


function  jp=jprime (n,x) 

%  Computes  the  first  derivative  of  the  Spherical  Bessel 
%  function. 

%  Written  by  Joseph  L.  Day 
jp=sqrt (pi/ (2*x) ) * (besselj (n- . 5,x) - 
(n+1) *bessel j (n+ . 5 , x) / (x) ) ; 


function  hp=hprime (n,x) 

%  Computes  the  first  derivative  of  the  Spherical  Bessel 
%  function. 

%  Written  by  Joseph  L.  Day 
hp=sqrt (pi/ (2*x) ) * (besselh(n- . 5, 2 , x) - 
(n+1) *besselh(n+.5,2,x) / (x) ) ; 


function  jn=sbessj (n,x) 

%  Computes  the  Spherical  Bessel  function 

%  Written  by  Joseph  L.  Day 
jn=sqrt (pi/ (2*x) ) *besselj (n+.5,x) ; 
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APPENDIX  B:  THREE -BODY  ADDITION  THEOREM  CALCULATION 

The  solution  for  the  mutual  radiation  impedance  of  the 
three-body  problem  as  described  in  the  Chapter  III  part  B 
using  the  Addition  Theorem  is  provided  in  this  appendix.  We 
start  by  writing  the  three  radial  velocities  as  expressed  in 
Equation  (2).  The  pressure  from  each  of  the  three  spheres 
(before  we  set  sphere  three  as  acoustically  hard)  can  be 
expressed  in  the  form  of  Equation  (3) . 

As  before  we  utilize  the  Addition  Theorem  to  write  the 
pressures  from  spheres  2  and  3  relative  to  sphere  one  (see 
Equation  (4)).  This  gives  us  the  expression  for  P. ^  and  P* . 
The  boundary  conditions  on  the  surface  of  each  of  the 
spheres  leads  to  the  following  three  equations. 


v,=— j-[p,+p'i+p,j] 

COpo  Or  rra 

v,=-^-[p.+p;+pI\ 

0)po  or  r2=a 

v,=-^-[pi+ Pl+ 

0)po  or  r2=a 

Then  we  obtain  three  big  equations  that  are  similar  to 
Equations  (7)  and  (8)  each  with  an  additional  term.  For 


example  looking  at  Equation  (7),  the  additional  term  will  be 
much  like  the  term  in  the  parenthesis  but  will  involve  the 
coefficients  for  the  third  sphere  (A3). 
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Following  a  similar  approach  as  used  in  Appendix  A,  we 
put  our  Equation  (7) -type  expressions  in  matrix  form.  We 
can  write  the  following  matrix  relation. 


7 

B\  B\ 

-  - 

cl 
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I  Bl 

-> 

A 
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where  A= 

A2 

Bl  I  _ 

_  _ 

Hi  .« 

to 

1 _ 

- - 1 

CO 

As  before  we  will  truncate  the  infinity  value  to  K  =  6 
gives  expression  of  the  following  sizes. 


- 

~A’s 

C1 

A’s~ 

-  ■  - 

C1 

147^147 

141x1 

— 

c2 

c3_ 

solving 

147x147 

/ 

c2 

c3_ 

We  have  solved  for  the  three  sphere  coefficients  and 
are  now  ready  to  complete  our  evaluation  of  the  mutual 
radiation  impedance.  Using  Equation  (13)  with  these  A- 
coefficients  I  can  solve  for  the  radiation  impedance.  Now  I. 
can  write  sphere  one's  total  radiation  impedance  as  Zrl  =  Zu 
+  Z12*V2/V1+Z13*Vi/V1.  Now  I  make  sphere  three  acoustically 

hard  by  setting  the  radial  velocity  amplitude  in  my  computer 
code  to  zero.  Next  I  solve  for  the  total  radiation 
impedance  by  setting  sphere  one  and  two's  radial  velocity 
amplitudes  equal  to  one.  I  store  this  result  and  then 
calculate  sphere  one's  self-radiation  impedance  by  setting 
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sphere  two's  radial  velocity  amplitude  to  zero.  I  subtract 
this  last  result  from  the  total  radiation  impedance,  and 
that  provides  my  mutual  radiation  impedance  for  sphere  one 
due  to  all  three  spheres.  The  plotted  data  is  calculated  by 
incrementing  the  range  of  sphere  three  as  it  moves  away  from 
the  two- sphere  axis. 

Here  is  the  MATLAB  computer  code  used  to  generate  these 
results.  All  supporting  functions  have  already  been 
provided. 


MATLAB  Program  for  3 -body  mutual  radiation  impedance 
%  Written  by  Joseph  L.  Day 

%  Mutual  impedance  using  the  Addition  Theorem 
%  Two  spheres  radiating,  third  (hard)  sphere  moving  away 
%  Last  updated  August  16,  1999 
format  long 

K=6;  %  Truncated  integer  vice  infinity 

%  A  useful  number 
%  waveNumber . 

%  for  all  three  spheres 
%  Angle  from  2 ' s  origin  to  1 ' s 
%  Defines  the  mode  for  each  sphere 

%  Monopole  position  for  matrices 
%  another  monopole; 

%  like  a  monopole 
%  spheres  velocity  amplitudes 
%  speed  of  sound  in  the  water 
%  density  of  the  water 

templ=eye (K2 ) ; temp2=zeros (K2 , K2 ) ; 

Big=[ tempi  temp2  temp2;temp2  tempi  temp2;temp2  temp2  tempi]; 

%  setup  some  space  for  future  calculation 

count=l;  %  initialize  a  counting 

process 

phicheck=0; 


K2=(K+1) A2; 
k=2 ; 

radius =0 . 5 ; ka=k* radius ; 

theta=0;phi=pi; 

nl=0;  ml=0;  n2=0;  m2=0; 

n3=0;  m3=0; 

positional; 

position2=l; 

position3=l; 

V1=1;V2=1;V3=0; 

speed=1490; 

rho=1000; 
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for  d=pi/2 : .25:24  %  Distance  with  0.25  meter 

increment 
dl2=pi; 

countla=l ; count2a=l ; count3a=l ; %  initialize  counts . 
countlb=l;count2b=l;count3b=l; 

Btemplda=zeros(l,K2) ;Btempldb=zeros(l,K2) ; 
Btemp2da=zeros(l,K2) ;Btemp2db=zeros (1,K2) ; 
Btemp3da=zeros(l,K2) ;Btemp3db=zeros(l,K2) ; 
for  t=0 : 6 

for  s=-t: t 

ptempla=abs (t-nl) ;ptemplb= 

abs (s-ml) ;ptemplc=max(ptempla,ptemplb) ; 
ptemp2a=abs (t-n2) ;ptemp2b= 

abs ( s-m2 ) ; ptemp2c=max (ptemp2a , ptemp2b) ; 
ptemp3a=abs (t-n3) ;ptemp3b= 

abs(s-m3)  ;ptemp3c=max(ptemp3a,ptemp3b)  ; 

count3=l ; count4=l ; count5=l ; count6=l ; count7=l ; count8=l ; 
for  pl=ptemplc : 2 : t+nl  %for  B(lto2) 

M=s-ml  ; 
if  pl<0 

P=  (-1) ''abs  (M)  *fact  (pl-abs  (M) )  / 

fact (pl+abs (M) ) *legendre (pi, cos (theta) ) ; 

else 

P=legendre (pi , cos ( theta) ) ; 

end 

Btempla=a(nl,pl,t,ml,s)*shank2(pl,dl2); 
Btemplb=P(l+abs (M) ) *exp (i*M*phi) ; 

Btemplc ( 1 , counts ) =Btempla*Btemplb; 
count3=count3+l ; 

end 

Btemplda (1, count la )=sum( Btemplc ) *jprime(nl,ka) / 

hprime (nl , ka) ; 

countla=countla+l ; 
phil3=pi/2+atan(pi/2/(d) ) ; 
thetal3=3*pi/2 ; 

for  pl=ptemplc : 2 : t+nl  %for  B(lto3) 

M=s-ml; 
if  pl<0 

P= (-1) Aabs (M) *fact (pl-abs (M) ) / 

fact (pl+abs (M) ) *legendre (pi, cos (theta!3) ) 

else 

P=legendre (pi , cos ( thetal3 ) ) ; 

end 

dl3=sgrt ( (pi/2) A2+ (d) A2) ; 

Btempla=a (nl ,pl, t,ml, s) *shank2 (pi, dl3 ) ; 
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Btemplb=P (1+abs (M) ) *exp (i*M*phil3 ) ; 

Btemplc (1, count4) =Btempla*Btemplb; 
count4=count4+l ; 

end 

Btempldb  (1,  countlb)  =sum (Btemplc)  *jprime  (nl,ka)  / 

hprime (nl , ka ) ; 

countlb=countlb+l ; 

for  p2=ptemp2c : 2 : t+n2  %for  B(2tol) 

M=s-m2; 
if  p2<0 

P= (-1) Aabs  (M)  *fact  (p2-bs  (M) )  /fact  (p2+abs  (M) )  * 
legendre (p2 , cos (theta -pi) ) ; 

else 

P=legendre (p2 , cos ( theta-pi )) ; 

end 

Btemp2a=a (n2 ,p2 ,  t,m2 ,  s)  *shank2  (p2 ,  dl2 )  ; 

Btemp2b=P (1+abs (M) ) *exp (i*M* (phi+pi) ) ; 

Btemp2c  ( 1 ,  count5 )  =Bteinp2a*Btemp2b; 
count5=count5+l ; 

end 

Btemp2da (1,  count2a)  =sum(Btemp2c)  *jprime (n2 ,ka)  / 

hprime  (n2 ,  ka )  ; 

count2a=count2a+l ; 
phi23=pi/2-atan (pi/2/ (d) ) ; 
theta23=3*pi/2 ; 

for  p3=ptemp3c :2 : t+n3  %for  B(2,3) 

M=s-m3; 
if  p3<0 

P= (-1) Aabs (M) *fact (p3-abs (M) ) / 

fact (p3+abs (M) ) *legendre (p3 , cos ( theta23 ) ) ; 

else 

P=legendre (p3 , cos ( the ta2 3 ) ) ; 

end 

d23=sqrt ( (pi/2) A2+ (d) A2) ; 

Btemp3a=a (n3 ,p3 , t ,m3 , s) *shank2 (p3 , d23 ) ; 
Btemp3b=P( 1+abs (M) ) *exp(i*M* (phi23) ) ; 

Btemp3c (1, count 6) =Btemp3a*Btemp3b; 
count6=count6+l ; 

end 

Btemp2db(l,  count2b)  =s\am(Btemp3c)  *  jprime  (n3 ,  ka)  / 

hprime  (n3 , ka)  ; 

count2b=count2b+l ; 

for  p3=ptemp3c : 2 : t+n3  %for  B(3,l) 

M=s-m3 ; 
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if  p3<0 

P=  (-1)  Aabs  (M)  *fact  (p3-abs  (M)  )  / 
fact (p3+abs (M) ) *legendre (p3 , cos ( thetal3-pi) ) 
else 

P=legendre (p3 , cos (thetal3-pi) ) ; 

end 

d23=sqrt((pi/2)A2+(d)A2); 

Btemp3a=a(n3,p3/  t,m3,s)*shank2  (p3,dl3)  ; 
Btemp3b=P ( 1+abs  (M) )  *exp  (i*M*  (phil3+pi) )  ; 

Btemp3c  ( 1 ,  count7 )  =Btemp3a*Btemp3b; 
count7=count7+l ; 

end 

Btemp3da (1, count3a) =sum(Btemp3c) * j prime (n3 ,ka) / 

hprime (n3 , ka) ; 
count3a=count3a+l; 
for  p3=ptemp3c:2 : t+n3  %for  B(3,2) 

M=s-m3 ; 
if  p3<0 

P= (-1) Aabs (M) *fact (p3-abs (M) ) / 
fact (p3+abs (M) ) *legendre (p3 , cos (theta23-pi) ) 
else 

P=legendre (p3 , cos ( theta23-pi) ) ; 

end 

Btemp3a=a (n3 ,p3, t ,m3 , s) *shank2 (p3 , d23 ) ; 
Btemp3b=P( 1+abs (M) ) *exp (i*M* (phi23+pi) ) ; 

Btemp3c (1, count8) =Btemp3a*Btemp3b; 
count8=count8+l ; 
end  ^ 

Btemp3db (1, count 3b) =sum(Btemp3c) * jprime (n3 , ka) / 

hprime (n3 , ka) 

count3b=count3b+l ; 

end 

phicheck=phicheck+ . 25 ; 
end 

Big (posit ionl , 2*K2+1 : 3*K2 ) =Btempldb;  %  Placement 
Big(positionl,K2+l:2*K2)=Btemplda; 

Big (K2+position2 , 1 : K2 ) =Btemp2da ;  %  of  vectors. 

Big (K2+position2 , 2*K2+1 : 3*K2 ) =Btemp2db; 
Big(2*K2+position3, 1 :K2) =Btemp3da; 
Big(2*K2+position3,K2+l:2*K2)=Btemp3db; 

C (3*K2 , 1) = zeros; 

C (positionl, 1) =-i *rho*speed*Vl /hprime (nl/ka) ; 
C(position2+K2,l)=-i*rho*speed*V2/hprime(n2,ka); 

C (position3+2*K2 , 1) =-i*rho*speed*V3 /hprime (n3 , ka) ; 
A=Big\C; 


Temporary=sum(conj (A(K2+1:2*K2,1) ) ' . *Btemplda*sbessj (nl,ka) ) 
*hprime(nl,ka) /jprime(nl,ka)  ; 

Temporary2=sum(conj (A(2*K2+1 : 3*K2 , 1) ) ' . *Btempldb*sbessj (nl,k 
a) ) *hprime(nl,ka) /jprime (nl,ka) ; 

Zrl=4*pi*radiusA2/Vl* (A(l, 1) *shank2 (nl,ka) +Temporary+Tempora 
ry2 )  ; 

Zrvec tor (count) =Zrl; 
count=count+l ; 

end 

Zrlscaled=Zrvector/ (4*pi*radiusA2*rho*speed) ; 
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APPENDIX  C:  T-MATRIX  CALCULATION 

Upon  obtaining  the  scattered  pressures  from  the  ATILA 
code  for  the  user  defined  mesh,  I  then  calculated  the 
Transition  Matrix  (T-Matrix)  for  the  user-defined 
transducer.  For  this  analysis  we  used  a  thin  spherical 
shell  for  our  transducer. 

Before  I  spell  out  my  method  for  calculating  the  T- 
matrix,  I  should  provide  the  reader  with  some  background 
concerning  this  T-matrix  (Ref.  3) .  First  suppose  the 
incident  pressure  to  an  object  are  standing  spherical  waves 
as  represented  by 

p(r,0,(t>)=  jn(kr)Qmn{$,<p). 

We  use  this  T-matrix  to  transform  from  our  incident 
pressure  to  the  scattered  pressure  from  the  object  subjected 
to  this  incident  pressure.  Basically,  the  T-Matrix 
describes  the  scattering  property  of  each  unique  scatter  (or 
in  this  case  transducer)  by  using  a  discrete  basis  set  of 
spherical  harmonics.  I  determine  the  T-Matrix  for  my 
spherical  shell  by  analyzing  the  scattering  characteristics 
of  this  shell  subject  to  a  "standing"  wave.  .  Mathematically 
speaking,  the  T-matrix  transforms  our  incident  standing 
waves  to  scattered  outgoing  waves  in  the  form  of  Hankel 
functions  shown  below. 
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p  (r, e, (f)  =  hn  [kr)Qmn  {d,  <p). 

The  following  equation  indicates  the  utility  of  the  T- 
Matrix.  The  Rn  values  represent  the  reflection  coefficients 
from  the  spherical  shell  and  constitute  elements  of  the  T- 
Matrix. 

T-Matrix  •  v  ' 

yZlPjp(kr)QmP(0^)  =>  ^Rnhn(kr)Q 

n  p  .  .  n  , 

" - - - V - '  ■  .  ' - — - - ' 

S  tan  ding  Wave  Radiating  Wave 

Here  is  how  I  used  these  ideas  to  develop  my  T-matrix 
calculation.  Starting  with  the  output  data  from  ATILA  which 
is  basically  Itf-number  of  data  points  in  the  following  form: 
for  i  =  1  to  N  we  have  xir  y.,  z.,  and  R.  Then  for  all  N- 

points  I  need  to  calculate  the  following  product 

h, (r,)£2l K<p, ) = h. ir, xk  )* P'A^e, y x £*. . 

I  will  do  this  product  up  to  the  quadrupole,  that  is 
nine  times  for  each  n  and  m  combination  (n=0  and  m=  0,  n=l 
and  m=-l,  ...  n=2  and  m=2)  .  So  for  example  if  we  have  100 

data  points,  the  number  of  calculated  complex  numbers  is  900 
(100  data  points  up  to  the  quadrupole).  The  resulting 
matrix  is  of  the  form: 
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calculated  complex  numbers 

'  t  ' 

A 

(i unknowns ) 

'  T  ' 

Pressure 

(fromATILA) 

= 

100  by  9 

(9x1) 

i 

(100x1) 

i 

We  solve  for  the  unknown  A-coefficients  in  MATLAB  by 
the  following  equation  [A]  =  [calculated  numbers ] \ [ P]  .  The 
forward  divide  sign  in  MATLAB  is  a  solution  for  the  A  vector 
in  a  least  squares  sense  to  this  over-determined  system  of 
equations  [calculated  numbers]* [A]  =  [P] .  The  effective 
rank  of  the  calculated  numbers  matrix  is  determined  from  the 
QR  decomposition  with  rank  pivoting.  The  relationship  of 
the  T-matrix  with  the  incident  pressure  amplitude  vector  [B] 
and  the  A-coef ficients  is  as  follows 


Boo 

-Aoo 

■“  “ 

Box 

Ao-i 

T  -  matrix 

Bio 

= 

Ao 

Bn 

•  ^ 

_ i 
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Since  we  obtain  bur  results  from  ATILA  for  each  n  and  m 
combination  we  solve  for  each  column  of  the  T-matrix  in  the 
following  way.  For  example,  we  obtain  ATILA' s  first 
scattering  pressure  results  when  we  model  the  incident 
pressure  with  the  standing .waves  using  n= 0  and  m= 0.  We  then 
use  the  scattering  pressures  to  solve  for  the  A- 
coef f icients .  Then  we  solve  for  the  first  column  of  the  T- 
matrix  by  solving  for  the  following  equation. 


V 

Am 

0 

Ao-i 

0 

= 

Ao 

0 

_  Al.22  _ 

Proceeding  this  way  we  can  solve  for  all  9-columns  of 
the  T-matrix  by  changing  the  n  and  m  values  of  the  incident 
pressure  accordingly.  One  final  comment  about  this 
procedure  to  calculate  the  T-Matrix  concerns  the  form  of  the 
incident  wave.  Because  of  certain  limitations  with  the 
ATILA  code  we  used  a  "traveling"  wave  vice  the  "standing" 
wave.  The  form  of  the  incident  pressure  wave  possible  is  a 
Hankel  function  instead  of  the  spherical  Bessel  function. 
Because  of  this,  a  correction  factor  to  the  diagonal 
elements  of  0 . 5* (diagonal  element  -  1)  is  necessary. 

The  following  text  provides  the  MATLAB  program  codes 
that  I  used  to  calculate  the  T-matrix.  Also  provided  are 
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the  user  defined  functions  used  in  support  of  the  main 
program  and  the  analytical  calculation  for  the  T-matrix  for 
thin  spherical  shells. 


Main  MATLAB  Program : 

%  Joseph  L.  Day 
%  Calculation  of  the  T-matrix 
%  last  updated  July  21,  1999 

%  Part  1:  Import  ATILA  data  and  calculate  the  A  coefficients 
format  long 

data=getdata ( ' temp5 . dat ' ) ; 

%  temp5.dat  is  a  text  file,  tab  deliminated  with  the  first 
%  line  as:  number  of  rows,  one-space,  and  number  of  columns. 
%  The  data  as  follows:  column  1:  x;  col.  2:  y;  col.  3:  z; 

%  col.  4:  real (pressure) ;  and 
%  column  5:  Imaginary (pressure) . 
x=data ( : , 1 ) ; y=data ( : , 2 ) ; z=data ( : ,  3 )  ; 
p=data(:,4)+i.*data(:,5); 

k=2*pi*474/1492 . 94;  %  since  ka=l  and  a=.5meters 

count=l;  %  initializes  my  count, 
for  N=0 : 2 

for  M=-N:N 

Ptemp2=zeros (1, length (x) ) ; 
for  n=l : length (x) 

.  r=sqrt ( (x(n) ) A2+ (y (n) ) A2+ (z (n) ) A2) ; 
theta=acos (z (n) /r) ; 
phi=atan2 (y (n) , x (n) ) ; 
if  M<0 

P= (-1) ^abs (M) *fact (N-abs (M) ) / 

fact (N+abs (M) ) *legendre (N, cos (theta) ) ; 

else 

P=legendre (N, cos (theta) ) ; 

end 

Ptempl=shank2 (N,k*r) *P(l+abs(M) ) *exp (i*M*phi) ; 
Ptemp2 (n) =Ptempl; 

end 

Ptemp2= (conj (Ptemp2) ) ' ;  %  make  column  vector 
eval ( [ ' Ptemp3 '  num2str (count)  ' =Ptemp2 ; ' ] ) ; 
count=count+l ; 

end 

end 
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cal=[Ptemp31  Ptemp32  Ptemp33  Ptemp34  Ptemp35  Ptemp36  Ptemp37 
Ptemp38  Ptemp39] ; 

%  cal  is  a  100by9  matrix  that  will  be  used  to  determine 

%  the  A  coefficients 

A=cal\p; 

%  This  final  "A"  value  provides  one  column  of  the  T-matrix. 

%  At  this  point  I  would  apply  the  correction  factor  since 
%  we  used  a  traveling  wave  instead  of  the  standing  wave. 

%  The  n  and  m  value  determines  which  column  of  the  T-matrix. 


Supporting  user  defined  functions : 
function  A=getdata (name) 

%Function  to  read  free-fo.rmated  data  A  (matrix  or  vector) 
%  Written  by  Charles  W.  Therrien 
if  all (name  -='.') 

name= [name, ' . dat 1 ] ; 
end 

[ ft, message] =f open (name, ' r' ) ; 
if  ft<0 

error (message) 

end 

[bf ,N] =f scanf (ft, ' %g' ) ; 
fclose (ft) ; 
if  bf (1) ==N-1 

A=zeros ( 1 , bf ( 1 ) ) ; 

A ( :  )  =bf  (2  :N)  ; 
elseif  bf (1) *bf (2) ==N-2 
A=zeros (bf (2) ,bf (1) ) ; 

A(  : ) =bf (3 :N) ;  ' 

A=A.  ’  ; 
else 

A=bf . 

fprintf ( t 1  Data  length  does  not  match  count . \n ’,. .. 

’  File  contents  is  being  returned  as  a  vector.An']) 

end 


function  hn=shank2 (n, x) 

%  Computes  the  Spherical  Hankel  function  of  the  second  kind. 

%  Written  by  Joseph  L.  Day 

hn=sqrt (pi / ( 2 *x) ) *besselh (n+ . 5 , 2 , x) ; 
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function  y=fact(n) 

%  Computes  the  factorial  of  n 
%  Written  by  Joseph  L.  Day 
y=prod(l :n) ; 


Program  to  calculate  the  Analytical  T-matrix  for  a  thin 
spherical  shell ; 

%  Written  by:  Joseph  L.  Day 

%  Reference:  ONR  Report  -  AY98  Thesis  (Ref  5) 

%  last  updated  April  6,  1999 

%  Computation  of  the  T-matrix  diagonal  elements  using 
%  the  analytically  derived  formulas  found  (Ref  5,  page  5) . 
format  long 

%  is  the  radius  of  the  sphere  in  meters 
%  is  the  shell  thickness 
%  is  the  bulk  modulus  in  Pascals 
%  is  the  Poisson’s  ratio  shell's  material 
%  is  the  angular  frequency  (474  Hz) 

%  is  the  density  of  the  solid  [kg/mA3] 

%  is  the  density  of  the  fluid  [kg/mA3] 

%  is  the  sound  speed  of  the  fluid  [m/s] 

%  is  the  wave  number,  in  this  case  ka=l 
cp=sqrt (E/ (rhos* (l-vA2 ) ) ) ;  %  plate  wave  speed 
omega =w*a/cp; 

M=1 ; 

for  N=0 : 2 
n=N; 

for  m=-n:n 

b=omegaA2-v-n* (n+1) +1 ; 
c=omegaA2-2* (1+v) ; 
d=n* (n+1) * (1+v) A2 ; 

Z(M)=i*h*cp*rhos/(a*omega)*((b*c-d)/b); 
jn=sqrt(pi/(2*k*a))*besselj(n+.5,k*a); 
hn=sqrt(pi/(2*k*a))*besselh(n+.5,2,k*a); 
jhat=sqrt (pi/ (2*k*a) ) * (besselj (n- . 5, k*a) 

- (n+1) *besselj (n+ . 5 , k*a) / (k*a) ) ; 
hhat=sqrt (pi/ (2*k*a) ) * (besselh (n- . 5, 2 ,k*a) 

- (n+1) *besselh(n+.5,2,k*a) / (k*a) ) ;  : 
numerator=i*Z (M) * jhat+rhof *cf * jn; 


a=0 . 5 ; 
h=0 . 01 ; 
E=215*10A9 ; 
v=0 . 33 ; 
w=2  *pi*474 ; 
rhos=7500 ; 
rhof=1000 ; 
cf =1492. 94; 
k=w/cf ; 
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denom=i*Z (M) *hhat+rhof *cf *hn; 

T  (M)  =-numerator/denom; 

M=M+1; 

end 

end 

results [ (1:9) '  real(T)'  imag(T) '  abs(T)'  (angle (T) *180/pi) ' ] 
%  result  is  the  output:  column  1 :  index  number  (n) ; 

%  column  2:  real (T-matrix) ; 

%  column  3:  image (T-matrix) ;  column  4:  amplitude; 

%  and  column  5:  phase  (degrees) . 
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